Digital Quantum Simulation with Rydberg Atoms 
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We discuss in detail the implementation of an open-system quantum simulator with Rydberg states 
of neutral atoms held in an optical lattice. Our scheme allows one to realize both coherent as well 
as dissipative dynamics of complex spin models involving many-body interactions and constraints. 
The central building block of the simulation scheme is constituted by a mesoscopic Rydberg gate 
that permits the entanglement of several atoms in an efficient, robust and quick protocol. In 
addition, optical pumping on ancillary atoms provides the dissipative ingredient for engineering 
the coupling between the system and a tailored environment. As an illustration, we discuss how 
the simulator enables the simulation of coherent evolution of quantum spin models such as the two- 
dimensional Heisenberg model and Kitaev's toric code, which involves four-body spin interactions. 
We moreover show that in principle also the simulation of lattice fermions can be achieved. As an 
example for controlled dissipative dynamics, we discuss ground state cooling of frustration-free spin 
Hamiltonians. 

PACS numbers: 03.67.-a, 05.30.Rt, 76.30.Mi, 32.80.Ee 



I. INTRODUCTION 

Simulating the evolution of many-body quantum sys- 
tems on classical computers is a complex task, which is 
believed to be intrinsically beyond the capabilities of clas- 
sical computation [1]. While in the classical case the 
number of degrees of freedom scales linearly with the 
particle number the computational complexity in treat- 
ing interacting many-particle quantum systems increases 
drastically due to an exponential growth of the Hilbert 
space dimension. One way to overcome this difficulty is 
to mimic the behavior of the quantum system of interest 
in an analogue physical system whose degrees of freedom 
are well accessible and controllable [H, Q . Initializing 
this system in a desired quantum state and measuring 
its properties after a given time is then effectively equiv- 
alent to having performed a simulation of the quantum 
evolution. Such quantum simulators are currently devel- 
oped for several physical platforms (see f° r a recent 
overview) , ranging from atomic systems [5|-[8j, trapped 
ions (914121 , im plementations based on nuclear magnetic 
resonance [l3|, Q > to photonic devices [HI, [l6| . 

One paradigmatic system for the simulation of many- 
body quantum physics are gases of ultracold alkali atoms 
trapped in optical lattices, where the underlying Hamil- 
tonian parameters such as the hopping rate or onsite in- 
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teraction strength energies of atoms can be controlled 
and tuned externally (see Fig. [1]) [I?). Since the first 
observation of a Mott-insulator to superfluid quantum 
phase transition enormous progress has been made in 
controlling these systems [18], including, very recently, 
the demonstration of single-site addressability [HI, [20| . 

Analogue quantum simulators are often purpose built, 
i.e., they are particularly well-suited for simulating spe- 
cific classes of many-body quantum systems. While the 
afore-mentioned ultracold atoms in optical lattices typ- 
ically realize bosonic or fermionic Hubbard-type Hamil- 
tonians with short-range interactions, trapped ions for 
instance naturally offer the possibility to study interact- 
ing spin models. In general, it is difficult to use analogue 
quantum simulators for the study of many-body systems 
with interactions, that differ considerably from the nat- 
ural, physically present (one- and two-body) interactions 
underlying the quantum simulator. For the simulation of 
exotic models involving certain constraints and higher- 
order interactions, these terms are usually created per- 
turbatively. This is typically associated by fine tuning 
problems, small effective energy scales and hence slow 
dynamics [2lj . 

These problems can be overcome by switching to a cir- 
cuit based model in a digital, which is sketched in Fig. [TJd. 
Here, the state of the system is encoded in qubits. The 
Hamiltonian time evolution is effectively created by a se- 
quences of quantum gates which act on these qubits, i.e. 
concatenating certain elementary gates will amount to 
the action of an effective time-evolution operator whose 
structure can be tailored with great flexibility. It has 
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FIG. 1: Analogue vs. digital quantum simulator, a: In an 
analogue simulator the interactions of the simulated model are 
typically closely related to the physical interactions that un- 
derlie the system which implements the simulator. In the case 
of cold atoms these are for example optical lattice potentials 
or short-ranged interatomic interactions, which can be used 
to simulate Bose and Fermi-Hubbard type Hamiltonians. b: 
In the digital case the Hamiltonian evolution is implemented 
by sequences of quantum gates acting on arrays of qubits, 
which are e.g. encoded internal states of atoms. The simu- 
lated Hamiltonian H e R can be vastly different from the inter- 
actions governing the underlying physical system. It allows 
for example the implementation of exotic many-body Hamil- 
tonians involving (higher-order) n-body interactions, which 
are very different to those normally encountered in ultracold 
atomic systems. 



been shown that, provided a universal set of quantum 
gates is available, such digital quantum simulator can be 
used to simulate the dynamics of any many-body Hamil- 
tonian with short-range interactions efficiently [3|- Fur- 
thermore it is possible to also efficiently simulate general 
(Markovian) dissipative dynamics, by including dissipa- 
tive reset operations on ancillary qubits [221-24] . Hereby, 
the engineering of a controlled coupling of the system to 
an artificially tailored environment offers the possibility 
of realizing dynamics for the dissipative preparation of 
entangled states and quantum phases [25], and cl osely 
related, quantum computation based on dissipation [24{ . 
Recently, we have developed a physical implementation 
of an open-system quantum simulator based on neutral 
atoms arranged in an optical lattice [26]. Alternative 
methods for the simulation of spin systems have recently 
been discussed: the ground state cooling for the toric 
code and a non-abelian topological phase using a single 
control atom moving through the lattice and interacting 
via two-qubit gates with the system spins has been pro- 
posed [27j . This method requires after an initial measure- 
ment procedure a second correction step, which removes 
the entropy from the system. In addition, the simula- 
tion of the coherent time evolution and the preparation 
of thermal states for the toric code using a stroboscopic 
method has been discussed by Herdman et al [28]. On 
the experimental side minimal instances of spin plaque- 
tte models have recently been implemented with trapped 
ions [l2| and photons [29|, [30| . 

In this work we review and extend such simulation ar- 
chitecture combining the coherent time evolution as well 



as dissipative terms for interacting spin systems, and in 
addition provide novel results on the simulation of lattice 
fermions. The central building block of our simulator is 
a mesoscopic Rydberg gate which relies on electromag- 
netically induced transparency (EIT) and the strong and 
long-ranged interaction of neutral atoms excited to Ryd- 
berg states [3l|. We furthermore show that by including 
optical pumping on ancillary atoms as a dissipative ingre- 
dient enables the implementation of open-system many- 
body dynamics. This dissipative dynamics can be used to 
perform efficient ground state cooling for a large class of 
spin models. In general one can expect that the cooling 
of any frustration-free Hamiltonian can be achieved. 



II. SIMULATION OF COHERENT DYNAMICS 

A. Setup and general scheme 

In the specific setup we have in mind ultracold atoms 
are trapped in a deep optical lattice in a Mott-insulator 
state with a single atom per site. These atoms are used 
to encode the qubits of our digital simulator in differ- 
ent electronic ground states. The spacing between the 
lattice sites can be on the order of up to a few mi- 
crometers. Experimentally, lattices which grant single 
site laser addressability have been demonstrated in Refs. 

As sketched in FigJTJ} the temporal evolution in a 
digital quantum simulator is achieved by concatenating 
quantum gates. In practice such gates are realized by 
carefully timed laser pulses and their interplay with state 
dependent interactions of the trapped atoms. The gen- 
eral aim is to simulate Hamiltonians of the form 



H 



(1) 



Here the hk are quasi- local Hamiltonians that govern the 
interaction of degrees of freedom located in the vicinity 
of the k-th lattice site. An appropriate sequence of gates 
will thus implement the time-evolution operator 



U (r) = exp 



(2) 



Note that the simulation time r is in general different 
from the real-time t. While the idea is simple the prac- 
tical implementation of such a scheme bears difficulties. 

Owed to the finite range of physical interactions, gate 
operations that can be efficiently implemented are usu- 
ally quasi- local. The time-evolution operator U(r), how- 
ever, contains highly non-local terms which in case of 
non-commuting hk cannot simply be decomposed into 
products of local ones. In practise one therefore tries 
to approximate U (r) by a sequence of quasi-local gates. 
Such approximation scheme is provided by the Suzuki- 
Trotter decomposition [34|. Here one approximately im- 
plements the global time-evolution operator over a time 
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step of length r by decomposing it into a product of time 
evolution operators, which correspond to the individual, 
quasi-local terms hk of the Hamiltonian: 



U(t) « l[exp[-i(T/h)h k ]. 



(3) 



Non-commutativity of terms hk leads to errors, which 
are bounded and can be controlled by the length of the 
time-step and reduced by choosing more sophisticated 
expansion schemes that are a generalization of eq. (j3j) 
(see, e.g., Ref. Q). 

While in theory any universal set of quantum gates 
allows the efficient approximation of arbitrary unitaries 
and therefore the dynamics of arbitrary Hamiltonians 
with short-range interactions, the experimental realiz- 
ability strongly demands for an implementation of the 
time-evolution with few and robust gates. While single- 
qubit gates are usually straight-forwardly implemented 
the challenge lies on creating two-qubit or even many- 
qubit gates as depicted in Fig. [TJd. This key requirement 
is met by a mesoscopic Rydberg gate which is detailed in 
the next subsection. It allows the implementation of an 
entangling multi-qubit gate on a microsecond timescale 
and with only three laser pulses, independently on the 
number of atoms involved in the multi-qubit gate. It 
therefore promises the flexibility, speed and robustness 
that is necessary for an efficient digital implementation 
of quantum spin models with exotic interactions. 



B. The mesoscopic Rydberg gate 

The purpose of this section is to review the basic prop- 
erties of the mesoscopic Rydberg gate presented in [31] . 
which constitutes the central building block of our en- 
visioned digital quantum simulator. The mechanism un- 
derlying the gate operation makes use of a two-photon in- 
terference phenomenon known as Electromagnetically In- 
duced Transparency (EIT) [35] . The general setup for the 
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FIG. 2: Setup for the mesoscopic quantum gate. A single 
control atom can be addressed independently of N ensemble 
atoms. Laser excitations induce a Rydberg interaction be- 
tween control and ensemble atoms, leading to the realization 
of a mesoscopic quantum gate. 

implementation of the quantum gate is shown in Fig. [2j 
We consider a single control atom and N ensemble atoms. 



For our setup we assume single-site addressability as it 
has recently been demonstrated by several experimental 
groups [HI, HO, [H, [33(| . The logical (qubit) states of the 
control atom are two hyperfine ground states denoted by 
|0) and |1). The logical states of the ensemble atoms are 
named \A) and \B). In spite of the different labeling it is 
in practice not necessary to distinguish between the con- 
trol and the ensemble atoms - an example for this will be 
given in Sec. Ill C II 

The mesoscopic Rydberg gate uses state-dependent in- 
teractions between Rydberg atoms [36l-[38| to realize a 
Controlled-NOT^ (CNOT^) gate, which is defined by 



G=|0)(0| c 



N 



N 



|l)(l|c.® 



(4) 



where - depending on the state of the control qubit - the 
state of all N target qubits is left unchanged or flipped. 
Here, af\A) i = \B) i and of IB), = [A),. 

To illustrate the underlying mechanism we introduce 
additional internal levels in both control and ensemble 
atoms which will be used to physically implement the 
gate operation (see Fig. [3]): The control atom has an 
auxiliary Rydberg state |r) that can be coupled to the 
hyperfine state |1) by a laser. In the ensemble atoms we 
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FIG. 3: Atomic level structure and external laser couplings for 
the mesoscopic gate. The states |1) and \r) in the control atom 
are coupled by a laser with Rabi frequency Q r - The weak laser 
fields Q p (t) drives a Raman transitions from \A) to \B) in the 
ensemble atoms, (a) For the control atom in |0) the Raman 
lasers and the strong coupling laser (Rabi frequency Q c ), cou- 
pling \P) to the \R) state, are in two-photon resonance, (b) 
For the control atom in |r) the Rydberg interaction shifts the 
\R) level away from the two-photon resonance. 

employ two additional levels. First, there is a coupling 
characterized by a time-dependent Rabi frequency Q p (t) 
between the hyperfine ground states \A) and \B) and 
an intermediate non- Rydberg \P) level, which still has 
a low principle quantum number such that interactions 
with the \r) level of the control atom are negligible. Fur- 
thermore, we make use of a Rydberg state \R) in each 
ensemble atom that is coupled to the intermediate \P) 
state with a laser of Rabi frequency Q c . The external 
laser fields are chosen such that there is a large detuning 
A from the \P) level, such that this state is only virtually 
populated. However, the hyperfine ground states and the 
Rydberg state are in two-photon resonance. 

With this setup the gate operation is performed by a 
sequence of three laser pulses that is depicted in Fig. 2J 



0: 

FIG. 4: Laser pulse sequence for the mesoscopic gate consist- 
ing of an initial tt pulse on the control atom, an adiabatic 
Raman transfer in the ensemble atoms, and a second tt pulse 
on the control atom. 

We start by applying a tt pulse on the control atom which 
transforms its qubit state a\0) + /3\1) to a\0) + i/3\r). In 
the next step we perform a conditional adiabatic Raman 
transfer in the ensemble atoms from \A) to \B) via the 
intermediate \P) state. To this end we apply a smoothly 
varying pulse profile £l p (t), which is chosen such that it 
realizes an effective 7r-pulse between \A) and \B). Finally, 
a second 7r-pulse is applied to the control atom. 

In the following we study the consequences of this pulse 
sequence on the ensemble atoms in the two cases in which 
the control atom is in |0) or |r), respectively. The full dy- 
namics then follows by taking the superposition accord- 
ing to the coefficients a and /3. Let us for simplicity first 
assume that the ensemble atoms do not interact with 
each other; consequences of non- vanishing interactions 
will be discussed later. Then, the dynamics of the en- 
semble atoms reduces to the product of the independent 
evolution of a single ensemble atom. For a large detun- 
ing A we may adiabatically eliminate the \P) level and 
obtain the effective Hamiltonian 

no 2 

tfeff = ^ [x 2 |+X+l + (1 + V)\R\R\ + x (\+XR\ + h.c.)] 

(5) 

Here, |+) = (\A) + \B))/y/2 is the symmetric superposi- 
tion of the two hyperfine ground states and x = V2Q P /Q C 
defines the relative strength of the probe laser Q p to the 
coupling laser Q c . Note that during the second laser 
pulse (Raman transfer) Q c and therefore also x are time- 
dependent functions. The interaction term V is in fact 
state-dependent and accounts for the state of the con- 
trol atom: In an ideal situation we have V = for the 
control atom in |0) (see Fig. [3^) while for the control 
atom in |r) (see Fig. [3b) we have a dominant Ryd- 
berg interaction, i.e., V = oo. The antisymmetric state 
|— ) = (\A) — \B))/y/2 is a zero energy eigenstate of the 
Hamiltonian. This dark state will thus be unaffected by 
the dynamics. 

Let us now look at the situation in which the control 
atom is in state |0) (Fig. [3^) and all ensemble atoms 
are in \A): Here the first and the last 7r-pulse shown in 
Fig. [4] have no effect. We have V = and we find that 
Hamiltonian (|5)) possesses in addition to | — ) a second 
zero energy dark state, 

W = (l + X i )- 1 /*[\+)-x\R)], (6) 

which for t = corresponds to the |+) state. The only 



non-zero eigenstate has the energy E2 = (l + x 2 ) 
and is thus energetically separated from the dark state 
manifold. During the Raman pulse (see Fig. |4]) the sys- 
tem will adiabatically follow the zero energy dark states 
for weak coupling lasers with x <C 1 and smooth laser 
pulse shapes Q p (t). Thus it will follow the dark state 
\d) = (1/V2) [|-) + which starts and ends in \A). 
The ensemble is hence effectively transparent for the Ra- 
man laser. Imperfections of this adiabatic passage arise 
from Landau-Zener transitions to the non-zero energy 
eigenstate [31[. 

In the case of the control atom starting in 1 1) the first 
7r-pulse will effectuate a transfer to the Rydberg state |r) 
and the strong interaction between the Rydberg levels V 
will change the outcome of the Raman laser sequence. 
For the sake of simplicity we assume V = 00. Here the 
Rydberg levels of the ensemble atoms will not take part 
in the dynamics as they are far off-resonant (see Fig. [3)3). 
Here the time evolution of the ensemble atoms follows the 
Hamiltonian H = M2 2 /(4A) x 2 |+X+|. Then, by choosing 
the pulse shape of Vt p {t) (x = x{t)) such that f x 2 (t)dt = 
7r the system will undergo the transformation (or Raman 
transfer) 

|->-H-> . 1+) ->-!+>■ (7) 

Expressing this transformation in the original states \A) 
and I B) results in 

\A)^-\B) , \B)^-\A), (8) 

which is the desired operation up to a trivial phase factor, 
which can be corrected by choosing suitable phases of the 
laser fields. Finally, the second 7r-pulse de-excites the 
control atom. 

Combining the two scenarios outlined above estab- 
lishes a way to control a NOT operation on the qubit 
states of the ensemble atoms conditioned on the state of 
the control atom, effectively realizing the CNOT^ gate 
(jlj) . In the next subsections we will make extensive use of 
this gate when implementing digital quantum simulation 
schemes for spin models (with many-body interactions). 

Before we proceed, however, we want to briefly dis- 
cuss issues relevant to the experimental implementation. 
First of all, we have neglected in our considerations the 
interaction among ensemble atoms. This is in general un- 
justified since the ensemble atoms are in the course of the 
gate sequence excited to Rydberg states which strongly 
interact. In a situation in which the control atom is ini- 
tially in 1 1) and hence is excited to a Rydberg state this 
does not constitute a major problem since the Rydberg 
state of the ensemble atoms is shifted far out of resonance 
and therefore is not excited. In the opposite situation 
(control atom in state |0)) the ensemble-ensemble inter- 
action is however expected to modify the working of the 
gate considerably. Indeed, one finds that in this case the 
state of the ensemble atoms can no longer be described 
by a tensor product of dark states [31]. Instead one finds 
that the initial state, e.g. \A N ) is written as a sum of 



5 



dark and "grey" states which acquire a dynamical phase 
shift (with respect to the dark states) that in the limit 
of infinite ensemble-ensemble interaction is proportional 
to Nmdix[x 2 (t)] [3l|| . In order to keep the fidelity of the 
gate high this shift has to be kept small. This means the 
higher the number of qubits that are to be entangled, 
the smaller the ratio x = \/2ft p /ft c . However, at the 
same time the Raman transfer has to be carried out at a 
time much shorter than the lifetime of the atomic Ryd- 
berg states, that is typically on the order of 50 /is. This 
requires a very strong coupling and also strong Raman 
lasers. One can show that the gate as presented here 
can implement an entangling operation on a timescale of 
~ 1 /is. In Ref. [3l| we have explicitly shown that fideli- 
ties > 99% can be achieved for N = 3 ensemble atoms 
located at a distance ~ 2/im from the control atom. The 
main limitation to the fidelity is given by the available 
laser power. In order to achieve high fidelities also for 
a larger number of ensemble atoms one has to ensure 
y/Nx <C 1. This choice suppresses deteriorating effects 
caused by the interaction among ensemble atoms excited 
to Rydberg states [39|, Ho|. Further consequence that 
arises from gate imperfections on the desired quantum 
simulation are analyzed towards the end of Sec. Ill C II 




FIG. 5: Lattice model for Kitaev's toric code, consisting of 
two sublattices involving plaquette operators A p and star 
terms B s . 



Fig. [5]). Besides its initially envisioned use as a quan- 
tum memory, this model has recently received consid- 
erable attention in the context of quantum simulation 
[27l I29I Hi| . The operators A p and B s are stabilizer op- 
erators with eigenvalues ±1. The model can be solved 
exactly, as all plaquette and star terms of the Hamilto- 
nian mutually commute. The global ground state is 
at the same time the ground state of each of the operators 
A p and B s : 



C. Digital simulation of spins and fermions 

In the following, we will use the mesoscopic Rydberg 
gate as building block for the digital quantum simulation 
[26j of spin Hamiltonians. We will at first discuss Kitaev's 
toric code whose Hamiltonian contains many-body inter- 
action. Despite the seemingly complicated structure this 
model is rather simple to implement as its Hamiltonian 
contains no non-commuting parts. The second example 
concerns the two-dimensional Heisenberg model and as a 
third example we discuss the digital simulation of lattice 
fermions. The latter can be mapped onto an effective 
spin Hamiltonian with six-body interaction terms. 



1. Kitaev's toric code 

In Kitaev's toric code model, spins are located on the 
links of a two-dimensional square lattice and interact via 
four-body interactions [4l[. This model is paradigmatic 
for a whole class of so-called stabilizer codes [HI, |43| . Its 
Hamiltonian is given by 



B.W) = IV) 



(10) 



(9) 



with the "plaquette terms" A p ^ = cr x a x a x a x being the 

product of four Pauli spin matrices and "star" terms 
(i) 

Bs = cr z a z a z a z defined in an analogous manner, see 
Fig. [5j In this notation the index j labels the plaque- 
ttes/stars and the four a-operators act on the spins lo- 
cated at the corners of the corresponding square (see 



for all plaquettes and stars, respectively. For periodic 
boundary conditions on a torus the stabilizers satisfy the 
relations Yl A p = 1 and Yl s B s = 1. For a system of 
N atoms there are N — 2 independent stabilizers. Con- 
sequently, the ground state manifold of the system will 
be four- fold degenerate. For an experimentally more ac- 
cessible situation of flat two-dimensional lattice struc- 
tures, the number of holes in the 2D lattice determines 
the ground state degeneracy. 

Excitations of the toric code Hamiltonian can be of two 
types: violations of the stabilizer constraints of either 
A p operators ("magnetic charges") or B p terms ("elec- 
tric charges" ) . They have an energy gap of AEq as every 
violation will affect two stars or plaquettes, respectively. 
In the following we will illustrate these excitations for the 
magnetic charges, but due to the symmetry of the Hamil- 
tonian the situation is identical for the electric charges. 

Flipping a single spin will create two magnetic charges 
located on adjacent plaquettes, see Fig. [6] By flipping 
a different spin on one of the adjacent plaquettes the 
excitations are effectively moved. The excitations are 
no longer quasi-local, but must be described by a string 
operator involving the path along which the charge has 
been moved. By flipping several spins we also may move 
a magnetic charge around an electric charge; due to the 
non-commutativity of <j x and <j z the state will eventually 
pick up a phase of tt. This behavior shows that the quasi- 
particles describing magnetic or electric charges neither 
have bosonic nor fermionic character as in both cases one 
would expect to recover the identity once the particle had 
been returned to its initial position. Hence, one calls such 
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particles with exotic statistics "anyons" related to their 
potential to pick up "any" phase under particle exchange. 
While the realizations of systems exhibiting anyonic ex- 
citations is interesting in itself, the anyon dynamics has 
direct consequences on the ground state cooling discussed 
in Sec. HITAl 




FIG. 6: Excitations in the toric code, (a) Flipping the spin 
indicated by the arrow will create two magnetic charges on 
the adjacent plaquettes. (b) Charges can be moved around by 
flipping further spins on the plaquettes containing the charges. 
The string operator characterizing the non-local excitation is 
shown as a solid line, (c) Moving a magnetic charge around 
an electric charge. 

Since the stabilizer operators A p and B s mutually com- 
mute no Trotter errors occur and the time-evolution op- 
erator for a time step r is exactly given by 

U = exp(—iHr/h) = exp(iEoA p r/h) ex.-p(iEoB s r/h), 

ps 

(11) 

and we need to focus only on the case of a single pla- 
quette or star. The extension to the entire lattice then 
follows naturally from iterating over all plaquettes and 
stars. Note that it is possible to parallelize many of these 
operations by partitioning the system into a few sublat- 
tices of piecewise independent atoms. 

During each timestep of the digital simulation we want 
the dynamics of the plaquette to be governed by the time 
evolution operator U p = ex.p(iEoA p t/h). In the regime 
of single site addressability, selecting a single plaquette is 
achieved by focusing the laser pulses which are required 
for the gate only on the atoms participating in the dy- 
namics of this plaquette. In Ref. [26j a simulation scheme 
was suggested where additional ancillary "control" atoms 
were used to effectively mediate four-body interactions 
between the four plaquette spins. Here, we propose an al- 
ternative approach which works without ancillary atoms. 
Instead one of the four plaquette atoms (as shown in Fig. 
[5j) will take the role of the control atom c. This situation 
requires a focussing of the Raman and control lasers ex- 
clusively on the remaining three atoms of the plaquette. 
This can either be achieved by an appropriate shaping 
of these laser beams or alternatively by addressing the 
three atoms sequentially between the excitation and de- 
excitation of the control atom. We can then decompose 
Up as 

U P = GU?(4>)G, (12) 
where G describes the mesoscopic Rydberg gate, 

G=|0)(0| c ® li + |l)<l| c ® of (13) 



and the single qubit rotation U^{(j>) = exp(icVr;f) imprints 
a phase of <fi = E$t/h onto the control atom. The sim- 
ulation of the star terms B s follows analogously by per- 
forming global 7r/2 rotation that interchange o x and <j z . 
More intuitively, the gate sequence that creates the effec- 
tive time-evolution under the Hamiltonian —EA P works 
as follows: The essential step is to realize that the inter- 
action 'sees' only whether a plaquette state has a positive 
or negative eigenvalue with respect to the operator A p . 
The actual configuration of the spins is not important. 
In the first step of the gate sequence (fT2j) the informa- 
tion of the state of three plaquette atoms is mapped on 
the control atom. Now the control atom is manipulated, 
i.e. a state-dependent rotation is carried out and the ini- 
tial mapping is reversed. The timescale of this process 
is essentially set by the speed of the mesoscopic Rydberg 
gate as single qubit rotations can be carried out quickly 
and reliably. Using experimentally realistic parameters, 
we find an effective interaction strength that can be on 
the order of several hundred kHz (26[. 

Let us now briefly investigate the effects of imperfec- 
tions on the quantum simulator. As mentioned above, a 
possible error source is the residual Rydberg interaction 
between the ensemble atoms during the mesoscopic gate, 
see Ref. (3l| for a detailed discussion. In the presence of 
such a coherent error we can write the mesoscopic gate 
as 

G'=|0)<0| c ®e i *« + |l)<l| c ® of. (14) 

Here, we have factored out the dependence on the im- 
printed phase cj) to allow for a consistent expansion. The 
perfect gate G then follows in the limit Q — >• 0, with the 
Hermitian operator Q acting on the entire plaquette ex- 
cept for the control atom. Expanding U' p = G f U*((f))G f 
up to first order in <p we obtain 

U' p = l + 2i0Q|O)(O| c + i</>Ap + O(0 2 ), (15) 

which corresponds to the coherent dynamics under the 
modified Hamiltonian H' = iJ + (l c — (T*)Q. Conse- 
quently, imperfections can be naturally integrated into 
the framework of this digital simulator with incoherent 
errors leading essentially to a finite temperature for the 
simulated system jiBf . 

2. Heisenberg model 

Let us now discuss a scheme for the simulation of coher- 
ent dynamics according to the Heisenberg Hamiltonian 

H = -\ E 0w°? + <VM + J *°*°i) + h J2 a i 

i,3 i 

(16) 

Here, J x (J y ,J z ) denotes the coupling strength of x-type 
(y— , z— type) spin-spin interactions between neighbour- 
ing spins, and h denotes the strength of a single particle 
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term playing the role of an effective magnetic field in z- 
direction acting on all spins. In contrast to the toric code 
Hamiltonian discussed above, here not all terms in the 
Hamiltonian commute. Thus the coherent time evolu- 
tion has to be realized in a Trotter expansion with small 
time steps to keep Trotter errors from non-commuting 
terms small. 

Simulation of the dynamics due to the magnetic field 
term for a small time step r is straightforward as 
exp(— ihr J2 i of) can be implemented by a global rota- 
tion of all spins, realized e.g. by a corresponding AC 
Stark shift applied to all atoms. The pairwise interac- 
tion terms can be built up from two-atom Rydberg gates 
G and single-qubit rotations as follows. A small time 
step ex.-p(iJ x Tcrf(jj /2) = ex.p(iOcrf CTj/2) of pairwise x- 
type interactions between neighbouring spins i and j can 
be written as 

exp(i0<7?o?/2) 
= exp(-Mrcrf /4) exp(i0<7? a* /2) exp(i7rof /2) 
= exp(-Mrcrf /4) exp(iOa? /2) 

exp(-i6»aj(l - ct?)/2) exp^Trcrf /4) 
= exp(-Z7rof /4) exp^of /2) 

[|0X0U ® lj + |lXl|i ® exp(-*6>(7?/2)] exp(i7r<jf/4) 

(17) 



Up to single-qubit rotations, this corresponds to an en- 
tangling two-qubit Rydberg gate (term in square brack- 
ets), where atom i takes the role of the control qubit, and 
atom j undergoes a conditional spin flip. For 6 = ir this 
spin flip takes place with unit probability and the opera- 
tion reduces to the gate @ discussed above. To minimize 
Trotter errors in the simulation, small values <C 1 are 
required. This can be readily achieved by adjusting the 
pulse length of the second pulse in the pulse sequence 
of the gate. For shorter times and / or smaller inten- 
sities of this pulse, the Raman transfer between logical 
states \A) ■ and \B) ■ of the target atom corresponds to 
only a partial population transfer, and thereby directly 
realizes the entangling operation in the last line of (fT7|) . 
The implementation of y and z interaction terms in (fTSj) 
can be done by combining the described procedure with 
single-qubit rotations on both atoms. We note that for 
a fixed set of parameters in the Hamiltonian ([16]) it can 
be beneficial to seek shorter decompositions of the evo- 
lution operator on each pair of sites to reduce the simu- 
lation cost in terms of the required two-qubit entangling 
gates. Note that in contrast to previous proposals for the 
quantum simulation of the Heisenberg model with cold 
atoms [46|, H3 , the discussed implementation with Ryd- 
berg atoms can acheive energy scales up to J x ~ 100 KHz, 
which are limited by the available laser power and not by 
the microscopic coupling constants. 



3. Fermi- Hubbard model in two dimensions 

The Heisenberg model can also be seen as the limiting 
case of a more general Hamiltonian that also involves the 
motion of the particles, which is the fermionic version of 
the Hubbard model (48| . Currently, there is great interest 
in the quantum simulation of the Fermi-Hubbard model 
in two-dimensional systems because of the expected re- 
lations to high-temperature superconductivity. Most ap- 
proaches are centered around an analog simulation using 
ultracold fermionic atoms in optical lattices [7|,l8|, but the 
experimental requirements for reaching the temperature 
regime of magnetic ordering or even the regime of d-wave 
superfluidity remain very challenging. Here, we outline 
a different approach where lattice fermions are mapped 
on a spin Hamiltonian with many-body interactions that 
can be implemented using our digital quantum simulator 

& 

The Hamiltonian describing the single-band Fermi- 
Hubbard model for spin 1/2 particles on a square lattice 
is given by 

H = -t c L c jct + U 5^n n n u , (18) 

<ij>a i 

where c\ a creates a fermion at site i with spin a and 
^ia = c \ a c io- is the corresponding number operator. The 
parameter t describes a hopping of the fermions to ad- 
jacent sites, while U accounts for the interactions of two 
fermions on the same lattice site. The Heisenberg model 
discussed above follows in the limit U —> oo at half filling. 

Since we are dealing with fermionic particles, the dig- 
ital quantum simulator needs to incorporate fermionic 
statistics. So far, we have discussed how to create a uni- 
versal quantum simulator for spin interactions. While 
spin 1/2 particles have the correct fermionic anticommu- 
tator {cr~,cr z + } = li on-site, they have bosonic statistics 
[<r~, o-j~] = off-site. To overcome this, one has to apply 
a Jordan- Wigner transformation, which has the form 

(H = l( do-]o-- (19) 

4 = l ®o-*at (20) 

ch = \^-<>t)- (21) 

Here, we have to introduce an enumeration of the sites 
of the two-dimensional lattice. This can be done, e.g, by 
starting in the lower left corner, moving to the lower right 
corner, go up one site, move to the left, and so on until the 
entire lattice has been scanned over. Then, the presence 
of a fermion on site i corresponds to the presence of a 
down spin. Conversely, an empty site in the fermionic 
picture translates into a spin up particle. Note that for 
spin 1/2 fermions, one has to replace each fermion by two 
spin 1/2 particles to encode all four possible combinations 
on each site. 
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Let us now look at how the operators in the Fermi- 
Hubbard Hamiltonian transform. Clearly, the on-site in- 
teraction results in interactions of the form (1 — cr|L)(l — 
cr|, ). On the other hand, the hopping terms become 

c\crCi+kcT = ViaSi^i+kcr + ^a^i^i^ka ' (22) 

with the string operator S? k = cr^ acting on 

the sites between i and i + fc. If we consider hopping 
in the horizontal direction, we see that within our enu- 
meration scheme the string operators drop out and the 
resulting spin operators remain local. However, this is 
not true in the vertical direction. There, we are left with 
non- vanishing string operators that run across the entire 
width of the lattice and consequently describe a highly 
nonlocal interaction. 

This obstacle can be overcome by introducing auxiliary 
degrees of freedom as shown in Ref. [50|. Such degrees 
of freedom are constituted by a fermion field di that is 



prepared in the ground state of the Hamiltonian 

^aux = —V ^ P^/ 5 j/Pj/ + i^/_i, (23) 

with the mutually commuting projectors Pyj' = (cW + 
d\ /o .)(dj' a — dt) and {i f , j'} partitioning the lattice into 
directed graphs, see Fig. The ground state of H au ^ is 
given by the condition Py ^> = 1 for all i', j' . By replacing 
the vertical hoppings in the Fermi-Hubbard Hamiltonian 
by 

c\ a Ci+k* ^ Civd+kvPi't'+k (24) 

it is possible to turn them into local spin operators when 
the Jordan- Wigner transformation is applied. The price 
one has to pay for this decoupling is that the resulting 
Hamiltonian contains six-body interactions [50|. Includ- 
ing all degrees of freedom, the Hamiltonian is of the form 



i,3,° 



^ fx x y y \ / -i \ j + 1 x y / ^ \ ( ~\ % \ 



',3,o- 



i,3 



V a 2i,2j,a a 2i+l,2j+l,a°'2i' ,2j' ,a a 2i> '+l,2j> ,a a 2i> ' + 1,2/ ,a a 2i> '+l,2j> '+1,ct 
i,3,o- 

T 7" ^ -2 Z X X XX 

V a 2i+l,2j>l,a a 2i,2j>2,a a 2^2j / + l,cr a 2i / + l^ 
i,3,o- 

t/ \ ^ z z y y y y 

V 2-^ a 2i+l,2j,a°'2i+2,2j+l,(T Cr 2i / +l,2j' J a (J 2i' +2,2j' ,a (7 2i' +l,2j' + l J a (J 2i / +2,2j' + l,cr 



2i+2,2j>l,a a 2i / + l,2j / + l,cr a 2i / +2,2j / + l,a a 2i / + l,2j / +2,cr Cr 2i / +2,2j / +2,cr' 



(25) 



r 



where we have moved to a two-dimensional notation for — t{ci^<^ 2 ^ + h.c), which transform into 
the lattice. 

ft l,2 = ~t (^l,l,t a 2,l,T + a l,l,a a 2,l,t) °l',l',t 

, x x z 4- y y % 

= ~ * cr l,l,t cr 2,l,t Cr l / 5 l / ,t _ t<7 l,l,t Cr 2,l J t <J l / ) l / J t 

= fti + h 2 - (26) 



For an experimental realization one needs to imple 
ment four spin 1/2 particles on every site of the orig- 



The two terms hi and h<i can be implemented sequen- 
tially according to the Suzuki- Trotter formula Eq. (|3]). 
The time-evolution according to hi can be simulated by 



inal fermionic model. This could be realized either by the gate sequence Ui = U$ \ v A GU*{(j))GU$ \ v t , where 

stacking up the square lattice in four layers or by being U H is the Hadamard gate, which interchanges a x and 

able to address all particles per site individually, e.g., <j z and the phase shift (j) = tr/h relates the duration of 

by choosing different hyperfine states. Note that the each timestep r to the coupling constant t. The term hi 

six-body interactions can again be mapped on the fa- can be implemented analogously, by including local gates 

miliar \\ { <jf form by applying local tt/2 rotations on interchanging <j x and a y on the sites 1,1, t an d 2, 1,|- 



the spins involving of. As an example, let us demon- 



While the experimental requirements for the imple- 



strate the digital simulation of the Fermion hopping term mentations of non-commuting three-body and six-body 
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interaction terms is certainly challenging, this example 
illustrates the principle power of a digital quantum sim- 
ulator. By using specially tailored many-body gates it 
might be possible to group several terms of the Hamil- 
tonian together and thus reduce the number of required 
gate operations. 

III. COOLING INTO MANY-BODY GROUND 
STATES 

So far, we have discussed the coherent simulation of 
spin and fermion lattice models. What we have left aside 
is the preparation of the initial state for this simulated 
dynamics. Typically, one is not interested in dynamics 
of arbitrary initial states, but rather in the behavior of a 
certain class of states, often those with low energy with 
respect to the simulated Hamiltonian. One possibility 
would be to adiabatically follow the ground state from 
an experimentally accessible initial state; for example, in 
analog simulation of lattice models with cold atoms, the 
system is first cooled to quantum degeneracy, with the 
lattice being ramped up adiabatically afterwards. Here, 
we follow a different route: we engineer the interaction 
with a dissipative environment in such a way that the 
resulting dynamics cools the system into the many- body 
ground state of the Hamiltonian of interest [24-2§|, l5lf . 
With such a dissipative element the dynamics is no longer 
unitary, but can be described by a quantum master equa- 
tion for the system density operator p, which is of the 
form 

^P = ~\ [H,p] + ( c iP°l - \ {4^i p}) , (27) 

and where the rates ji control the strength of the dissi- 
pation. The goal is then to engineer the jump operators 




FIG. 7: Lattice model for the Fermi- Hubbard model for a sin- 
gle spin degree of freedom. Sites corresponding to the original 
fermions are colored in red, while the sites corresponding to 
the auxiliary fermion field are colored in blue. The arrows in- 
dicate the graphs along which the projectors Pij are defined. 




FIG. 8: Setup for the cooling of the toric code. Interstitial 
control atoms are shown in red, with their internal level struc- 
ture allowing optical pumping into the |0) state. 

Ci in such a way that the only stationary states of the 
master equation correspond to the groun dstates of the 
Hamiltonian of interest. Note that this Hamiltonian does 
not necessarily correspond to the Hamiltonian H gener- 
ating the coherent dynamics in the master equation. In 
fact, we will study the case of purely dissipative dynamics 
with H = in the following. 

A. Kitaev's toric code 

Let us focus again on the toric code whose coherent 
simulation was already discussed Sec. Ill C II In contrast 
to the implementation of the coherent dynamics we will 
no longer use one of the atoms of each plaquette or star 
as the control particle. Instead, we add control atoms to 
the interstitial spaces of our lattice as shown in Fig. [51 
In our cooling scheme these control atoms will be opti- 
cally pumped, and therefore provide the coupling to the 
dissipative environment needed for our quantum state 
preparation. 

As demonstrated in Sec. Ill C li the global ground state 
of the toric code is at the same time the ground state 
of each plaquette or star. Kitaev's toric code falls into 
the class of frustration-free Hamiltonians, where the en- 
ergy of each term can be minimized independently. We 
can therefore perform the cooling to the ground state 
of each plaquette or star. As in the coherent case we 
first focus on the case of a single plaquette A p . Here 
the local Hamiltonian is given by h p = —Eq A p with the 
spectrum of A p constituted by two eightfold degenerate 
sectors with eigenvalues +1 and —1. It is therefore con- 
venient to denote the states by |±1,A), where ±1 refers 
to the eigenvalue of A p and A labels the different states 
within the degenerate manifold. The preparation of the 
ground state sector can be achieved by pumping the pla- 
quette into any superposition or mixture of +1 eigen- 
states. This is realized by choosing a four-body quantum 
jump operator of the form 

c P = \°t (1 - A p ) , (28) 

where erf acts on an arbitrary spin i of the four plaquette 
spins. To understand the action of this jump operator, 
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it is instructive to split it into two parts. First, the "in- 
terrogation" part 1/2(1 — A p ) checks whether the system 
is already in the correct eigenstate. Applied to any +1 
eigenstate (A p = +1) the jump operator vanishes and 
the ground state manifold is left unchanged. However, 



for A n 



-1 the second "pump" part erf flips the sign 



of A p and consequently transforms any state | — 1,A) of 
the excited state manifold directly to the corresponding 
state 1 1 , A) in the ground state manifold. 

The implementation of this jump operator in terms 
of a gate sequence essentially follows this picture. The 
auxiliary control atom is initially prepared in the state 
|0) c . Then, the many-body eigenstate of A p = ±1 is 
mapped onto the control atom by the gate sequence 



S = RI{k/2)- 1 GRI{k/2), 



(29) 



where i?^(7r/2) = exp(— z7t<j^/4) is a local tt/2 rotation 
acting on the control atom and G is the mesoscopic Ryd- 
berg gate (jH). The mapping S can be described as 



|0)J+1,A> ^ |0)J+1,A> 
|0) C |-1,A) -> |1) C |-1,A) 



(30) 
(31) 



After this mapping we can therefore conditionally manip- 
ulate the many-body states of the plaquette by a condi- 
tional operation based on the state of the control atom. 
We perform a controlled spin flip onto one of the four 
system spins, given by 

U?(0) = |0X0| C 1 + |1X1| C exp(i0<jf ). (32) 

Here, the angle controls the probability with which a 
spin flip from the —1 to the +1 eigenspace is realized 
(see below). The two-qubit gate U?(9) can be imple- 
mented based on the mesocopic Rydberg gate. Finally, 
we reverse the mapping by applying the inverse gate se- 
quence S~ 1 (= S). Then, we find that the control atom 
remains in the |1) state every time Uf induces a spin flip, 
i.e. in general it remains entangled with the four plaque- 
tte spins. Consequently, before using the control atom for 
the next cooling step, optical pumping to the |0) state is 
required to reinitialize the control atom in |0) such that 
it factors out from the dynamics of the plaquette spins. 
It is this dissipative element that provides the necessary 
ingredient that allows one to remove entropy from the 
system. For <C 1 we can perform an expansion of the 
dynamical map describing the evolution of the density 
operator p, i.e., 



l(c pP cl-\{clc p ,p}^ +O(0 3 ) (33) 



dt P 



with the jump operators c p given in Eq. ([28]) and the 
cooling rate 7 = 6 2 /t. 

As demonstrated before, this discussion can be general- 
ized to the entire lattice system, with the jump operators 
for the site terms B s being given by c s = erf (1 — B s )/2. 
The cooling process can then be understood in the anyon 




FIG. 9: Numerical simulations of the dissipative state prepa- 
ration of the ground state of the toric code for N — 32 par- 
ticles. For times t large compared to the simulation timestep 
r, essentially all anyons are removed from the system and the 
ground state energy E — —NEo is reached asymptotically. 
By increasing the phase shift per timestep, the cooling effi- 
ciency can be enhanced. Each data point corresponds to an 
average over 1,000 realizations. 



picture as follows: Each spin flip incoherently moves an 
anyon to an adjacent plaquette or site, see Fig. [6j and 
whenever two anyons of the same type meet, they are 
annihilated and their energy is removed from the sys- 
tem. For larger values of the dynamics is given by a 
discrete version of the quantum master equation with the 
cooling being even more efficient. As can be seen from 
Fig. [9j the most efficient preparation of the ground state 
occurs for = tt. For this value of a dissipative move 
of an anyon to an adjacent plaquette takes place with 
unit probability. 



B. Frustration-free Hamiltonians 

The above analysis for ground state cooling of the toric 
code can be extended to a large class of interesting mod- 
els. In general, one can design jump operators, where 
any ground state of a frustration free Hamiltonian is the 
unique dark state. With a suitable choice, one would 
expect the cooling of any initial state into the ground 
state of the frustration free Hamiltonian in analogy to 
the toric code discussed above. An example for a spin 
liquid phase at the Roskhar-Kievelson point has been 
discussed in (26|. An important question for the exper- 
imental realization is the efficiency of such ground state 
cooling: within quantum information theory, the efficient 
cooling requires that the time evolution of an arbitrary 
initial state approaches the ground state exponentially 
with a characteristic cooling rate, which scales poly norm- 
ally in the system size. So far it has been demonstrated 
that an important subclass of frustration free Hamilto- 
nians, namely stabilizer states, can be cooled efficiently 
|24l l51j . However, it remains an open question whether 
the ground state of a general frustration free Hamiltonian 
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can be prepared efficiently with dissipative techniques. 
The toric code discussed above represents an example of 
a stabilizer state, which exhibits abelian topological or- 
der. However, a stabilizer formulation can also be derived 
for large classes of states exhibiting non-abelian topolog- 
ical order: an important example are the string net con- 
densates [52|. As a consequence, the dissipative ground 
state cooling illustrated for the toric code above allows 
one also to efficiently prepare ground states with highly 
non-abelian topological order. These can serve as the 
building block for a topological quantum computer. 

IV. CONCLUSION AND OUTLOOK 

In this work we aimed at discussing the implementa- 
tion of a digital quantum simulation architecture using 
Rydberg atoms in optical lattices. Our goal was further- 
more to outline schemes for the simulation of coherent 
and dissipative dynamics corresponding to (many-body) 
spin models. Recently, these concepts for the digital sim- 
ulation of open-system dynamics have been extended to 
systems of trapped ions [53]. In a remarkable experi- 
ment [l2|, a combination of single- and multi-qubit (en- 
tangling) gates and optical pumping has been used to 
simulate coherent four-body spin interactions and dis- 
sipative four-qubit stabilizer pumping, thereby demon- 
strating in a minimal system of one plaquette the elemen- 
tary building blocks required for future large-scale simu- 
lations of Kitaev's toric code and related models. While 
the underlying physical interactions of this trapped ions 
simulator naturally differ from the van-der-Waals inter- 
action between Rydberg atoms, the current experiments 
demonstrate the experimental feasibility of the digital ap- 
proach to quantum simulation of open-system dynamics 
in interacting many-body systems. Furthermore, they 



clearly illustrate the generic effect of how errors in the 
employed gates, as discussed in Sec. Ill C li affect the ac- 
tual simulated dynamics, and suggest that, ultimately, 
in future large-scale fault-tolerant quantum simulation 
quantum error correction techniques might have to be 
incorporated. 

For the Rydberg simulator architecture described in 
the present work, in principle all individual building 
blocks have been experimentally demonstrated individ- 
ually: The first entangling Rydberg gates have been re- 
cently realized in the laboratory for two atoms held in 
optical tweezers Ell- Important next steps will be 
to improve these Rydberg entangling operations and to 
incorporate them in larger optical trap arrays and / or 
optical lattices, where Mott insulator states involving 
many lattice sites can be prepared with high accuracy 
and single-site addressability has become available. In 
combination with the development of efficient multi-atom 
entangling Rydberg gates as the one reviewed here (3l| 
this appears to be a promising route towards large-scale 
quantum simulations of many-body spin models, which 
lie beyond the cope of classical simulation capabilities. 
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